void vscale(caja *bx)
{
	///Termostato empleando el "Momentum Scalling Method"

	///Variables externas: Ahi que ver de donde vienen y como generarlas
	double atemp = bx->temp;
	double temp=0.0;
	double taut=5000.0;
	////Variables de entrada del termostato
	double alpha;
	for(int j=0;j<3;j++)
	{
		for(int i=0;i<bx->npart;i++)
		{
			temp+= (bx->part[i].vel[j])*(bx->part[i].vel[j]);
		}
	}
	temp*=1.0/(3.0*bx->npart*kB);
	//alpha = sqrt(1.0+(atemp/temp-1.0)/taut);
	alpha = sqrt(atemp/temp);

	//alpha = sqrt(atemp/temp);
	//alpha=1.0;
	cout<<"alpha: "<<alpha<<endl;
	cout<<"atemp: "<<atemp<<"  "<<temp<<endl;
	bx->alpha= alpha;
}

void vscale2(caja *bx)
{
	///Termostato empleando el "Momentum Scalling Method"

	///Variables externas: Ahi que ver de donde vienen y como generarlas
	double atemp = bx->temp0;
	double temp=0.0;
	////Variables de entrada del termostato
	double alpha;
	for(int j=0;j<3;j++)
	{
		for(int i=0;i<bx->npart;i++)
		{
			temp+= (bx->part[i].vel[j])*(bx->part[i].vel[j]);
		}
	}
	//temp*=1.0/(3.0*bx->npart*kB);
	temp*=1.0/(3.0*bx->npart);
	alpha = sqrt(atemp/temp);
	for(int i=0;i<bx->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			bx->part[i].vel[j]*=alpha;
		}
	}
	//cout<<"atemp: "<<atemp<<"  "<<temp<<endl;
}

void dscale(caja *bx)
{
//	///Manostato modificando el volumen de simulacion y reescalando posicion de
//	///las particulas
//
//	///Variables externas: Ahi que ver de donde vienen y como generarlas
//	double beta, delta, apres,pres,taup, comp, facpot;
//	//////////////////
//      int nfi,nacc;
//      beta=pow(1.0-comp*delta*(apres-pres)/taup,1.0/3.0);
//
//      for(int i=0;i<3;i++)
//      {
//    	  bx->lbox[i]*=beta;
//    	  for(int j=0;j<bx->npart;j++)
//    	  {
//    		  bx->part[j].pos[i]*=beta;
//    	  }
//      }
//      bx->vol=bx->lbox[0]*bx->lbox[1]*bx->lbox[2];
}

void values(caja *bx, int nacc, bool printScreen)
{
	double epot, ekin, viri, pig, pint;

	//**************************Parametros a revisar****************************
	//double range;///alcance del potencial de interaccion
	//double facpot= 2.0*pi*range*range*range*eps*npart*pow(sigma,6)/(3.0*volume);
	double apres=1.0; //presion de entrada
	double facpot=0.0;///No se considera la correccion a la presion
	double pres,work,eges,rho,freeEn,freeEnH;
	bool kpres=false;//Activacion del manostato (testear)
	bool ktemp =false;//Activacion del termostato (escalado de velocidades) testear
	int lcrtl=1; //Frecuencia de correccion de temperatura y presion
	double temp;
	//**************************************************************************


	//Calculo de las energias
	epot = bx->epoten/bx->npart-facpot;
	ekin=bx->ekin/(bx->npart);
	///Control de la temperatura
	temp=2.0/(3.0)*ekin;
	///Virial para el Lennard-Jones (Sin correcciones)
	viri=0.5/bx->npart*bx->virial;
	///pressure according to Eq. 2.2.9 in Hansen and Mc Donalds
	///PIG (units of Kj/mol divide nm^3)
	pig=bx->npart*kB*temp/bx->vol;
	///VIRIAL expresado en Kj/mol (ri-rj*Fij sum i sum j>i)
	pint=-1.0/3.0*bx->virial/bx->vol;
	pres=pig+pint;

	//Imprimimos los resultados
	char name[50];
	sprintf (name, "%s", "resultados.txt");
	ofstream MiArchivo (name,ios::out|ios::app);
	if (MiArchivo.is_open())
	{
		MiArchivo<<nacc<<" "<<epot<<" "<<ekin<<" "<<epot+ekin<<" "<<temp<<" "<<viri<<" "<<pres<<" "<<pig<<" "<<pint<<endl;

		MiArchivo.close();
	}

	///Control del volumen
	work=1.0e-3*avogadro/bx->npart*apres*bx->vol;

	///Energia total por molecula en el paso nacc
	eges=ekin+epot;
	if(kpres)
		eges+=work;

	if(printScreen)
	{
		//Energia potencial y cinetica del sistema
		//cout<<"Epoten: "<<bx->epoten<<"  Ekin: "<<bx->ekin/kB<<endl;
		//Energia potencial y cinetica del sistema por particula
		cout<<"Epoten/N: "<<epot<<"  Ekin/N: "<<ekin<<endl;
		//Deriva del sistema
		cout<<"Deriva: "<<bx->vel[0]<<"  "<<bx->vel[1]<<"  "<<bx->vel[2]<<endl;
		//Virial
		cout<<"viri: "<<viri<<endl;
		//Presion
		cout<<"pres: "<<nacc<<"  "<<double(pres*bx->vol)<<endl;
		//Temperatura cinetica
		cout<<"Temperatura: "<<temp<<endl;
	}

	///densidad
	rho=bx->npart/avogadro/bx->vol;

	///Energia libre
	freeEn=eges;
	freeEnH=freeEn/(kB*temp);
	if(abs(freeEnH)>170.0)
		freeEnH=sign(170,freeEnH);

	///Control de la Temperatura y presion
	vscale2(bx);
	if(nacc%lcrtl==0)
	{
		if(kpres)
			dscale(bx);
		if(ktemp)
			vscale(bx);
	}
}
